module parameter_setting

    implicit none

    ! model parameters

    ! parameters given by a matlab code ###############################################################3

    double precision :: parameters(13)
    double precision :: parameters_ss(16)
    double precision :: parameters_integer(8)
    double precision :: parameters_solution(2)
    double precision :: init_guess(6)
    integer :: init_type

    double precision :: policy_parameters(4), tau_l, thre_l, thre_n, thre_pib

    double precision :: beta
    double precision :: nu
    double precision :: alpha
    double precision :: delta
    double precision :: lambda
    double precision :: chi0
    double precision :: chi1
    double precision :: rho_a
    double precision :: sigma_a
    double precision :: nd

    double precision :: psi
    double precision :: n0
    double precision :: gamma
    double precision :: sigma_rk

    double precision :: c_ss,y_ss,i_ss,k_ss,h_ss,n_ss,l_ss,rb_ss,rkhats_ss,a_ss,xi_ss,pr_ss

    ! welfare   
    double precision :: welfare_integer(2)
    double precision :: welfare_states(4)
    double precision :: utility, period_utility, expected_u, n1, k1, rb1, a1
    integer :: length,repeat,s
    double precision :: utility_compare

    ! simulation
    ! simulation
    double precision, allocatable :: av(:),kv(:),nv(:),rbv(:),xiv(:)
    double precision, allocatable :: cv(:,:),yv(:),iv(:),hv(:,:),lv(:),rhseev(:),pibv(:),kv_unreg(:),lv_unreg(:)
    double precision :: rkhat,rkhats,rkhat_next,rkhats_next
    double precision :: erk_next,pibv_next,erk_next_unreg
    double precision :: xi_next,c_next,h_next,y_next,rhsee_next,rv_next,a_next
    double precision :: z_temp
    double precision :: z_temp_log, cdf_log

    double precision, allocatable :: x0(:,:),x(:),x0_next(:,:),x_next(:),x0_reg(:,:),x_reg(:)
    double precision, allocatable :: rhsee_each(:),rhsee_derived(:)
    double precision :: dummy1,dummy2,rhsee_temp,cdf,pdf
    double precision :: mean_rhsee_norun,mean_rhsee_run
    double precision :: mean_rbar_norun ,mean_rbar_run
    double precision, allocatable :: rbar_derived(:)

    ! polynomials
    integer :: nstv,degp,npol ! number of state variables, degree of polynomial, size of polynomial matrix
    double precision, allocatable :: expmx(:,:),expmx_column(:),state(:)
    double precision, allocatable :: eta1_rhsee(:),eta2_rhsee(:),eta3_rhsee(:)
    double precision, allocatable :: eta1_rbar(:),eta2_rbar(:),eta3_rbar(:)

    ! trapezoid integration
    integer :: n_trapz
    double precision, allocatable :: z_trapz(:)

    ! other variables #################################################################################

    ! values for non-linear solver
    double precision, parameter :: tol_fun = 1d-8   
    double precision :: fval(1),temp(1)
    integer :: info

    ! policy function iteration and convergence
    integer :: npf
    double precision, allocatable :: shocks(:,:), shocks_column(:)
    integer :: t_max, iter_max, t_drop
    integer :: t, i, j, f, g, z, count

    ! to import data
    character(len=20) :: filename               

end module

! main program ##########################################################################################   

program baseline_model

    ! modules
    use parameter_setting
    use my_toolbox

    implicit none

    ! subroutines
    double precision, external:: residual_fn, residual_fn_unreg

    ! import parameters from text files #################################################################

    open (10, file='from_matlab_to_fortran/parameters.txt', status='old', action='read')
       read(10,*) parameters
    close(10)
    beta     = parameters(4)
    nu       = parameters(5)
    alpha    = parameters(6)
    delta    = parameters(7)
    lambda   = parameters(8)
    chi0     = parameters(9)
    chi1     = parameters(10)
    rho_a    = parameters(11)
    sigma_a  = parameters(12)
    nd       = parameters(13)

    open (10, file='from_matlab_to_fortran/ss_values.txt', status='old', action='read')
       read(10,*) parameters_ss
    close(10)
    c_ss      = parameters_ss(1)
    y_ss      = parameters_ss(2)
    i_ss      = parameters_ss(3)
    k_ss      = parameters_ss(4)
    h_ss      = parameters_ss(5)
    n_ss      = parameters_ss(6)
    l_ss      = parameters_ss(7)
    rb_ss     = parameters_ss(8)
    rkhats_ss = parameters_ss(9)
    a_ss      = parameters_ss(10)
    xi_ss     = parameters_ss(11)
    pr_ss     = parameters_ss(12)
    psi       = parameters_ss(13)
    n0        = parameters_ss(14)
    gamma     = parameters_ss(15)
    sigma_rk  = parameters_ss(16)

    open (11, file='from_matlab_to_fortran/parameters_integer.txt', status='old', action='read')
        read(11,*) parameters_integer
    close(11)
    npf      = nint(parameters_integer(1)) ! number of policy functions to be approximated
    nstv     = nint(parameters_integer(2)) ! number of state variables
    degp     = nint(parameters_integer(3)) ! degree of polynomial
    npol     = nint(parameters_integer(4)) ! size of polynomial matrix
    iter_max = nint(parameters_integer(5))
    n_trapz  = nint(parameters_integer(6)) ! number of grids for integration approximation
    t_max    = nint(parameters_integer(7))
    t_drop   = nint(parameters_integer(8))

    open (11, file='../no_policy_3rd/from_matlab_to_fortran/welfare_integer.txt', status='old', action='read')
        read(11,*) welfare_integer
    close(11)
    length = nint(welfare_integer(1))
    repeat = nint(welfare_integer(2))

    open (11, file='../no_policy_3rd/from_matlab_to_fortran/welfare_states.txt', status='old', action='read')
        read(11,*) welfare_states
    close(11)
    n1  = welfare_states(1)
    k1  = welfare_states(2)
    rb1 = welfare_states(3)
    a1  = welfare_states(4)

    allocate(expmx(nstv,npol))
    allocate(expmx_column(nstv*npol))
    open (11, file='from_matlab_to_fortran/polynomial_matrix.txt', status='old', action='read')
        read(11,*) expmx_column
    close(11)
    f=0
    do i=1,nstv
        do j=1,npol
            f=f+1
            expmx(i,j) = expmx_column(f)
        end do
    end do

    allocate(shocks(repeat,length),shocks_column(repeat*length))
!    open (11, file='C:\Users\hmats\Dropbox\Intel_Fortran\three_million_shocks_stdnormal.txt', status='old', action='read')
    open (11, file='../three_million_shocks_stdnormal.txt', status='old', action='read')
        read(11,*) shocks_column
    close(11)
    f=0
    do i=1,repeat
        do j=1,length
            f=f+1
            shocks(i,j) = shocks_column(f)
        end do
    end do

    allocate(eta1_rhsee(npol),eta2_rhsee(npol),eta3_rhsee(npol))
    allocate(eta1_rbar(npol),eta2_rbar(npol),eta3_rbar(npol))

    open (11, file='from_fortran_to_matlab/solved_eta1_rhsee.txt', status='old', action='read')
    read(11,*) eta1_rhsee
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta2_rhsee.txt', status='old', action='read')
    read(11,*) eta2_rhsee
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta3_rhsee.txt', status='old', action='read')
    read(11,*) eta3_rhsee
    close(11)

    open (11, file='from_fortran_to_matlab/solved_eta1_rbar.txt', status='old', action='read')
    read(11,*) eta1_rbar
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta2_rbar.txt', status='old', action='read')
    read(11,*) eta2_rbar
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta3_rbar.txt', status='old', action='read')
    read(11,*) eta3_rbar
    close(11)

    allocate(cv(repeat,length),hv(repeat,length))
    allocate(av(length),nv(length),kv(length),rbv(length),xiv(length))
    allocate(yv(length),iv(length),lv(length),rhseev(length),pibv(length),lv_unreg(length),kv_unreg(length))

    allocate(state(nstv))
    allocate(x0(nstv,npol),x(npol))

    open (11, file='from_matlab_to_fortran/policy_parameters.txt', status='old', action='read')
        read(11,*) policy_parameters
    close(11)
    tau_l    = policy_parameters(1)
    thre_l   = policy_parameters(2)
    thre_n   = policy_parameters(3)
    thre_pib = policy_parameters(4)

    ! repeat stochastic simulations ########################################################

    utility = 0d0

    do s = 1, repeat

    ! set initial states
    nv(1)  = n1
    kv(1)  = k1
    rbv(1) = rb1
    lv(1)  = k1/n1
    av(1)  = a1

    ! simulate for t_max periods
    do t = 2, length

        ! Variables determined independent of bank runs
        av(t)  = rho_a*av(t-1) + sigma_a*shocks(s,t-1)

        hv(s,t)  = ((1d0-alpha)/psi*exp(av(t))*kv(t-1)**(alpha))**(1d0/(alpha+1d0/nu))
        yv(t)  = exp(av(t)) * kv(t-1)**(alpha)*hv(s,t)**(1-alpha)

        rkhats = log ( rbv(t-1) * (1d0-1d0/lv(t-1)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
        rkhat  = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                + (1d0+nu)/(alpha*nu+1d0)*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t-1))
        xiv(t)  = exp(rkhat)/exp(rkhats)

        state(1) = log(kv(t-1))
        state(2) = log(rbv(t-1))
        state(3) = log(nv(t-1))
        state(4) = log(xiv(t))
        do i=1,nstv
            do j=1,npol
                x0(i,j) = state(i)**expmx(i,j)
            end do
        end do
        x  = x0(1,:)*x0(2,:)*x0(3,:)*x0(4,:)
        if ( xiv(t) >= 1d0 ) then
            pibv(t)  = ( exp(rkhat) + 1d0-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1d0)*nv(t-1) - nv(t-1)
            if ( pibv(t) > 0d0 ) then
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0
                rhseev(t) = exp( sum(x*eta1_rhsee) )
                rbv(t)    = exp( sum(x*eta1_rbar ) )
            else
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0
                rhseev(t) = exp( sum(x*eta3_rhsee) )
                rbv(t)    = exp( sum(x*eta3_rbar ) )
            end if
        else
            pibv(t)  = ( exp(rkhat) + 1d0-delta )*lv(t-1)*nv(t-1) &
                     - rbv(t-1)*(1d0+lambda*(1d0-gamma))*(lv(t-1)-1d0)*nv(t-1) - nv(t-1)
            nv(t) = n_ss*(1d0 - nd/100d0)
            rhseev(t) = exp( sum(x*eta2_rhsee) )
            rbv(t)    = exp( sum(x*eta2_rbar ) )
        end if

        cv(s,t) = rhseev(t)**(-1d0) + psi*hv(s,t)**(1d0+1d0/nu)/(1d0+1d0/nu)

        iv(t)  = yv(t) - cv(s,t)
        kv(t)  = (1d0-delta)*kv(t-1) + iv(t)

        lv(t)    = kv(t)/nv(t)
        erk_next = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                 + (1d0+nu)/(alpha*nu+1d0)*rho_a*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t))

        ! To obtain deposit rate R_bar, solve bank's FOC w.r.t. levarage
        temp(1) = rbv(t)
        call hybrd1 ( residual_fn, 1, temp, fval, tol_fun, info, t )
        rbv(t)  = temp(1)

        if ( info /= 1 ) then
            write(*,*) 'Equation not solved at period',t
        end if

!        if ( pibv(t) >= thre_pib .and. nv(t) > thre_n ) then
            
            lv_unreg(t) = lv(t)/(1d0-tau_l)
            kv_unreg(t) = lv_unreg(t)*nv(t)
            erk_next_unreg = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                           + (1d0+nu)/(alpha*nu+1d0)*rho_a*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv_unreg(t))
            temp(1) = rbv(t)
            call hybrd1 ( residual_fn_unreg, 1, temp, fval, tol_fun, info, t )
            rbv(t)   = temp(1)
            if ( info /= 1 ) then
                write(*,*) 'Equation not solved at period',t,', policy imposed.'
            end if

!        end if

        period_utility = beta**(t-2) * log ( cv(s,t) - psi*hv(s,t)**(1d0+1d0/nu)/(1d0+1d0/nu) )
        utility = utility + period_utility

    end do ! end of each period

    end do ! end of each simulation

    ! End of simulation #############################################################################

    expected_u = utility / repeat

    write(*,*) 'expected utility =',expected_u

    open (50, file='from_fortran_to_matlab/utility_compare.txt', status='replace', action='write')
        write(50,*) expected_u
    CLOSE(50)

end program

! residual function to obtain deposit rate rbv #########################################################
subroutine residual_fn(n_in, temp_in, residual, info_in, t_in)

    use parameter_setting
    use my_toolbox

    implicit none

    integer, intent(in) :: n_in, info_in, t_in
    double precision :: temp_in(n_in)
    double precision, intent(out) :: residual(n_in)
    double precision :: rbv_nls,rks_nls,z_nls,ex_nls,p_nls,pdf_nls,cdf_nls

    rbv_nls = temp_in(1)

    rkhats_next = log ( rbv_nls * (1d0-1d0/lv(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
    z_temp     = (rkhats_next - erk_next) / sigma_rk
    z_temp_log = (erk_next + sigma_rk**2d0 - rkhats_next) / sigma_rk
    call normp( z_temp    , cdf    , dummy2, pdf    )    ! obtain cdf and pdf
    call normp( z_temp_log, cdf_log, dummy1, dummy2 )    ! obtain cdf

    ! residual function to solve
    residual(1) = rbv_nls - &
                  ( exp(erk_next + sigma_rk**2d0/2d0) * cdf_log / (1d0-cdf) + (1d0-delta) &
                    - lambda * ( rbv_nls**2d0 * (1d0-gamma) * (1d0+lambda*(1d0-gamma)) * (lv(t)-1d0)/lv(t)**2d0 ) &
                             / ( rbv_nls * (1d0-1d0/lv(t))  * (1d0+lambda*(1d0-gamma)) - 1d0+delta ) &
                             /sigma_rk * pdf/(1d0-cdf) )

end subroutine
    
! residual function to obtain deposit rate rbv #########################################################
subroutine residual_fn_unreg(n_in, temp_in, residual, info_in, t_in)

    use parameter_setting
    use my_toolbox

    implicit none

    integer, intent(in) :: n_in, info_in, t_in
    double precision :: temp_in(n_in)
    double precision, intent(out) :: residual(n_in)
    double precision :: rbv_nls,rks_nls,z_nls,ex_nls,p_nls,pdf_nls,cdf_nls

    rbv_nls = temp_in(1)

    rkhats_next = log ( rbv_nls * (1d0-1d0/lv_unreg(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
    z_temp     = (rkhats_next - erk_next_unreg) / sigma_rk
    z_temp_log = (erk_next_unreg + sigma_rk**2d0 - rkhats_next) / sigma_rk
    call normp( z_temp    , cdf    , dummy2, pdf    )    ! obtain cdf and pdf
    call normp( z_temp_log, cdf_log, dummy1, dummy2 )    ! obtain cdf

    ! residual function to solve
    residual(1) = rbv_nls - &
                  ( exp(erk_next_unreg + sigma_rk**2d0/2d0) * cdf_log / (1d0-cdf) + (1d0-delta) &
                    - lambda * ( rbv_nls**2d0 * (1d0-gamma) * (1d0+lambda*(1d0-gamma)) * (lv_unreg(t)-1d0)/lv_unreg(t)**2d0 ) &
                             / ( rbv_nls * (1d0-1d0/lv_unreg(t))  * (1d0+lambda*(1d0-gamma)) - 1d0+delta ) &
                             /sigma_rk * pdf/(1d0-cdf) )
    
end subroutine
